Bernoulli numbers turn up in a surprising number of contexts. This page will document the most outstanding and make sense of why these numbers appear therein.

Contents

Definition of Bernoulli Numbers
Trigonometric Power Series
Sums of Products of Bernoulli Numbers
Sums of Powers of Integers
Bernoulli Polynomials
Euler-Maclaurin Summation
Relations to Zeta Functions

Definition of Bernoulli Numbers

These numbers are defined by the generating function

$\frac{t}{e^t-1}=\sum _{n=0}^{\infty }{B}_n\frac{t^n}{n!}$

The sum on the left-hand side is the inverse of an exponential series with powers shifted down by one. Rearranging the definition gives

$1=\sum _{k=0}^{\infty }{B}_k\frac{t^k}{k!}\times \sum _{m=1}^{\infty }\frac{t^{m-1}}{m!}=\sum _{k,l=0}^{\infty }{B}_k\frac{t^{k+l}}{k!(l+1)!}$

The first term on the right-hand side is constant, so that $B_0=1$ . The remaining powers must then have coefficients equal to zero. Picking out an arbitrary power gives

$\sum _{k=0}^n\frac{B_k}{k!(n-k+1)!}=0=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n+1}{k}\right)B_k$

This summation that equals zero is all important in detecting the presence of Bernoulli numbers, and will be used time and again in what follows.

Since binomial coefficients are positive by definition, some of the Bernoulli numbers must be negative for the total sum to equal zero. The numbers can be evaluated recursively using

$B_n=-\frac{1}{n+1}\sum _{k=0}^{n-1}\left(\genfrac{}{}{0ex}{}{n+1}{k}\right)B_k$

With the aid of a CAS or a bit of calculational persistence, the first few values are

$\begin{array}{l}{B}_1=-\frac{1}{2},B_2=\frac{1}{6},B_3=0,B_4=-\frac{1}{30},B_5=0,\\ B_6=\frac{1}{42},B_7=0,B_8=-\frac{1}{30},B_9=0,B_{10}=\frac{5}{66},\dots \end{array}$

All odd numbers after $B_1$ are identically zero. This is easy to show by separating the first two terms of the sum:

$\begin{array}{c}\sum _{n=0}^{\infty }{B}_n\frac{t^n}{n!}=1-\frac{t}{2}+\sum _{n=2}^{\infty }{B}_n\frac{t^n}{n!}=\frac{t}{e^t-1}\\ 1+\sum _{n=2}^{\infty }{B}_n\frac{t^n}{n!}=\frac{t}{e^t-1}+\frac{t}{2}=\frac{t}{2}\frac{e^t+1}{e^t-1}=\frac{t}{2}cot\frac{t}{2}\end{array}$

Since the final right-hand side is an even function, there can be no odd terms on the left-hand side. The coefficients of these terms are necessarily zero.

Even Bernoulli numbers alternate in sign, but this will not be important in what follows.

Trigonometric Power Series

Because of the generating function definition, Bernoulli numbers appear prominently in the power series of trigonometric functions other than sine and cosine. Consider first the hyperbolic cotangent:

$\begin{array}{l}cothx=\frac{e^x+e^{-x}}{e^x-e^{-x}}=\frac{e^{2x}+1}{e^{2x}-1}=1+\frac{2}{e^{2x}-1}\\ \phantom{cothx}=1+\frac{1}{x}\sum _{k=0}^{\infty }{B}_k\frac{(2x{)}^k}{k!}\\ cothx=\frac{1}{x}+\sum _{k=2}^{\infty }\frac{2^k{B}_k}{k!}{x}^{k-1}=\frac{1}{x}+\sum _{k=1}^{\infty }\frac{2^{2k}{B}_{2k}}{\left(2k\right)!}{x}^{2k-1}\end{array}$

The index of the sum can be rewritten in the last step because all the odd Bernoulli numbers in the sum are identically zero. The singular term can be combined with the sum but is kept separate for clarity.

For the circular cotangent, let $x\to ix$ and, noting that $sinhix=isinx$ , multiply by the i from the denominator of the left-hand side:

$cotx=\frac{1}{x}+\sum _{k=1}^{\infty }\frac{(-1{)}^k{2}^{2k}{B}_{2k}}{\left(2k\right)!}{x}^{2k-1}$

The hyperbolic tangent has a positive sign in its denominator, requiring a change of sign in the denominator of the the generating function. First note the partial fraction decomposition

$\frac{1}{e^{2t}-1}=\frac{1}{(e^t-1)(e^t+1)}=\frac{1}{2}[\frac{1}{e^t-1}-\frac{1}{e^t+1}]$

so that one has

$\frac{t}{e^t+1}=\frac{t}{e^t-1}-\frac{2t}{e^{2t}-1}=\sum _{k=0}^{\infty }{B}_k\frac{t^k-(2t{)}^k}{k!}$

The power series for the hyperbolic tangent is thus

$\begin{array}{l}tanhx=\frac{e^x-e^{-x}}{e^x+e^{-x}}=\frac{e^{2x}-1}{e^{2x}+1}=1-\frac{2}{e^{2x}+1}\\ \phantom{tanhx}=1-\frac{1}{x}\sum _{k=0}^{\infty }{B}_k\frac{(2x{)}^k-(4x{)}^k}{k!}\\ tanhx=\sum _{k=2}^{\infty }\frac{2^k(2^k-1)B_k}{k!}{x}^{k-1}=\sum _{k=1}^{\infty }\frac{2^{2k}(2^{2k}-1)B_{2k}}{\left(2k\right)!}{x}^{2k-1}\end{array}$

where again the index of the sum can be rewritten in the last step because all the odd Bernoulli numbers in the sum are identically zero.

For the circular tangent, let $x\to ix$ and divide by the i from the numerator of the left-hand side:

$tanx=\sum _{k=1}^{\infty }\frac{(-1{)}^{k+1}{2}^{2k}(2^{2k}-1)B_{2k}}{\left(2k\right)!}{x}^{2k-1}$

For the hyperbolic cosecant, use an intermediate partial fraction decomposition:

$\begin{array}{l}cschx=\frac{2}{e^x-e^{-x}}=\frac{2e^x}{e^{2x}-1}=\frac{1}{e^x-1}+\frac{1}{e^x+1}\\ \phantom{cschx}=\frac{1}{x}\sum _{k=0}^{\infty }{B}_k\frac{2x^k-(2x{)}^k}{k!}\\ cschx=\frac{1}{x}+\sum _{k=2}^{\infty }\frac{(2-2^k)B_k}{k!}{x}^{k-1}=\frac{1}{x}+\sum _{k=1}^{\infty }\frac{(2-2^{2k})B_{2k}}{\left(2k\right)!}{x}^{2k-1}\end{array}$

One more time the index of the sum can be rewritten in the last step because all the odd Bernoulli numbers in the sum are identically zero. The singular term is again kept separate for clarity.

For the circular cosecant, let $x\to ix$ and multiply by the i from the denominator of the left-hand side:

$cscx=\frac{1}{x}+\sum _{k=1}^{\infty }\frac{(-1{)}^k(2-2^{2k})B_{2k}}{\left(2k\right)!}{x}^{2k-1}$

Series for the hyperbolic and circular secants cannot be written in terms of Bernoulli numbers, but rather define a new set known as Euler numbers.

Sums of Products of Bernoulli Numbers

Consider the consistency of power series in the context of the two derivatives

$\frac{d}{dx}cothx=-{csch}^{2}x\frac{d}{dx}cschx=-cschxcothx$

Including singular terms in sums for conciseness, these two relations can be expressed as

$\begin{array}{l}\sum _{k=0}^{\infty }\frac{2^{2k}{B}_{2k}}{\left(2k\right)!}(2k-1)x^{2k-2}=-\sum _{k,l=0}^{\infty }\frac{(2-2^{2k})(2-2^{2l})B_{2k}{B}_{2l}}{\left(2k\right)!\left(2l\right)!}{x}^{2k+2l-2}\\ \\ \sum _{k=0}^{\infty }\frac{(2-2^{2k})B_{2k}}{\left(2k\right)!}(2k-1)x^{2k-2}=-\sum _{k,l=0}^{\infty }\frac{(2-2^{2k})2^{2l}{B}_{2k}{B}_{2l}}{\left(2k\right)!\left(2l\right)!}{x}^{2k+2l-2}\end{array}$

Picking out an arbitrary power of the independent variable on each side of the equation, these two expressions are equivalent to

$\begin{array}{l}(2n-1)\frac{2^{2n}{B}_{2n}}{\left(2n\right)!}=-\sum _{k=0}^n\frac{(2-2^{2k})(2-2^{2n-2k})B_{2k}{B}_{2n-2k}}{\left(2k\right)!(2n-2k)!}\\ \\ (2n-1)\frac{(2-2^{2n})B_{2n}}{\left(2n\right)!}=-\sum _{k=0}^n\frac{(2-2^{2k})2^{2n-2k}{B}_{2k}{B}_{2n-2k}}{\left(2k\right)!(2n-2k)!}\end{array}$

This will appear much simpler with the abbreviations

$S_1=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{2n}{2k}\right)B_{2k}{B}_{2n-2k}{S}_2=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{2n}{2k}\right)2^{2k}{B}_{2k}{B}_{2n-2k}$

Note that both sums exhibit an interchange symmetry between $k$ and $n-k$ . Keeping this in mind, the two equalities can be written

$\begin{array}{l}(1-2n)2^{2n}{B}_{2n}=(4+2^{2n})S_1-4S_2\\ (1-2n)(2-2^{2n})B_{2n}=2S_2-2^{2n}{S}_1\end{array}$

A bit of algebra establishes that $S_1=S_2$ , except for $n=1$ when the system is indeterminate. From either of these two relations, and for $n\ne 1$ , one then has

$\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{2n}{2k}\right)B_{2k}{B}_{2n-2k}=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{2n}{2k}\right)2^{2k}{B}_{2k}{B}_{2n-2k}=(1-2n)B_{2n}$

This is somewhat surprising: the extra power of two inside the second summation has no effect on the result. These summations can also be written without the first and last terms as

$\sum _{k=1}^{n-1}\left(\genfrac{}{}{0ex}{}{2n}{2k}\right)B_{2k}{B}_{2n-2k}=\sum _{k=1}^{n-1}\left(\genfrac{}{}{0ex}{}{2n}{2k}\right)2^{2k}{B}_{2k}{B}_{2n-2k}=-(2n+1)B_{2n}$

This result has been known since Euler, but perhaps this is a novel derivation.

Sums of Powers of Integers

Bernoulli numbers first arose in the context of summing powers of integers. To motivate the general form, first consider a trivial sum over zeroth powers:

$\sum _{k=1}^n{k}^0=\sum _{k=1}^{n}1=n$

The corresponding sum over first powers uses a trick of Gauss,

$\sum _{k=1}^{n}k=\frac{n}{2}(n+1)=\frac{1}{2}{n}^2+\frac{1}{2}n$

where the first and last terms are added and multiplied by half the number of such intermediate sums. This form remains an integer for n both even and odd.

The sum over squares of integers is most easily done using a trick of Knuth. Consider an identity for the sum of cubes:

$\begin{array}{l}\sum _{k=1}^n{k}^3+(n+1{)}^3=\sum _{k=0}^n(k+1{)}^3=\sum _{k=0}^n(k^3+3k^2+3k+1)\\ =\sum _{k=1}^n{k}^3+3\sum _{k=1}^n{k}^2+3\sum _{k=1}^{n}k+(n+1)\end{array}$

The sum over cubes cancels on both sides, but one can then solve for the sum of squares using previous results:

$\sum _{k=1}^n{k}^2=\frac{1}{3}[(n+1{)}^3-\frac{3}{2}n(n+1)-(n+1)]=\frac{1}{3}{n}^3+\frac{1}{2}{n}^2+\frac{1}{6}n$

This method can be used for all subsequent sums, as well as the previous two results.

These examples are enough to see that the final result begins with a power one more than that in the sum, descending to a single power. Thus assume that the sum over integers to an arbitrary power has the form

$\sum _{k=1}^n{k}^p=\sum _{k=1}^{p+1}{c}_k{n}^k$

For p = 0, this form is satisfied with c₁ = 1. Consider induction on the assumed form:

$\begin{array}{l}\sum _{k=1}^{n+1}{k}^p=(n+1{)}^p+\sum _{k=1}^n{k}^p\\ \phantom{\sum _{k=1}^{n+1}{k}^p}=(n+1{)}^p+\sum _{k=1}^{p+1}{c}_k{n}^k=\sum _{k=1}^{p+1}{c}_k(n+1{)}^k\end{array}$

Before expanding binomials, let $n\to n-1$ in the second line to simplify,

$n^p+\sum _{k=1}^{p+1}{c}_k(n-1{)}^k=\sum _{k=1}^{p+1}{c}_k{n}^k$

so that now the right-hand side can be easily eliminated:

$n^p+\sum _{k=1}^{p+1}\sum _{m=1}^k{c}_k(-1{)}^m\left(\genfrac{}{}{0ex}{}{k}{m}\right)n^{k-m}=0$

The indices k and m describe a triangular region in their coordinate plane. Diagonals in this triangular region correspond to the difference between the two indices. These differences run from zero, when the two indices are equal, up to the power p, when the indices differ the most. Working down in value from the uppermost diagonal, when the indices are equal, produces an new upper limit on the second index m. The double sum can thus be rearranged as

$n^p+\sum _{i=0}^p\sum _{m=1}^{p+1-i}{c}_{i+m}(-1{)}^m\left(\genfrac{}{}{0ex}{}{i+m}{m}\right)n^i=0$

One can now set coefficients of powers of n equal to zero for coefficient values. The highest power gives

$1-c_{p+1}\left(\genfrac{}{}{0ex}{}{p+1}{1}\right)=0\to c_{p+1}=\frac{1}{p+1}$

while the next highest power gives

$c_{p+1}\left(\genfrac{}{}{0ex}{}{p+1}{2}\right)-c_p\left(\genfrac{}{}{0ex}{}{p}{1}\right)=0\to c_p=\frac{1}{2}$

which is curiously independent of the power p. Both of these results are consistent with the examples at the beginning of this section. For the general case one has

$\sum _{m=1}^{p+1-i}{c}_{i+m}(-1{)}^m\left(\genfrac{}{}{0ex}{}{i+m}{m}\right)=0$

and the result just above suggest one should look at a descending index defined by $i+m=p+1-k$ . The general case then becomes

$\sum _{k=0}^{p-i}{c}_{p+1-k}(-1{)}^k\left(\genfrac{}{}{0ex}{}{p+1-k}{p+1-i-k}\right)=0$

where a possible overall minus sign is ignored as inconsequential. Now consider casting this in the form of the all-important summation for Bernoulli numbers:

$\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n+1}{k}\right)B_k=0$

One would think at first that the alternating nature of the sum in question would make this impossible, except for one fact: all odd Bernoulli numbers are zero except $B_1$ . Apart from explicitly adjusting the sign of the coefficient containing this value, one can ignore the alternating sign in the series.

One thus needs to have

$c_{p+1-k}\left(\genfrac{}{}{0ex}{}{p+1-k}{p+1-i-k}\right)\approx \left(\genfrac{}{}{0ex}{}{p-i+1}{k}\right)B_k$

Ignoring factors not dependent of the index of summation k, this gives

$c_{p+1-k}\approx \frac{B_k}{k!(p+1-k)!}\approx \left(\genfrac{}{}{0ex}{}{p+1}{k}\right)B_k$

Renaming the index and adjusting the normalization factor to match the first two values gives the result

$c_k=\frac{(-1{)}^{{\delta }_{k,p}}}{p+1}\left(\genfrac{}{}{0ex}{}{p+1}{k}\right)B_{p+1-k}$

The extra minus sign accounts for the one term in the alternating series that contains the nonzero odd Bernoulli number $B_1$ .

The final form of the sum over integers to an arbitrary power is thus

$\sum _{k=1}^n{k}^p=\sum _{k=1}^{p+1}\frac{(-1{)}^{{\delta }_{k,p}}}{p+1}\left(\genfrac{}{}{0ex}{}{p+1}{k}\right)B_{p+1-k}{n}^k$

and the leading terms are

$\sum _{k=1}^n{k}^p=\frac{1}{p+1}{n}^{p+1}+\frac{1}{2}{n}^p+\frac{p}{12}{n}^{p-1}-\frac{p(p-1)(p-2)}{6!}{n}^{p-3}+\dots$

where as noted $c_{p-2}\equiv 0$ , along every other subsequent coefficient, due to odd Bernoulli numbers being zero after the first. The coefficients here are naturally consistent with the examples at the beginning of this section. The linear term will not appear if multiplied by an odd Bernoulli number other than $B_1$ .

This form of the sum has a rather surprising implication,

$\begin{array}{c}\sum _{k=1}^n{k}^p=n^p+\sum _{k=1}^{n-1}{k}^p\\ \sum _{k=1}^{n-1}{k}^p=\sum _{k=1}^n{k}^p-n^p=\sum _{k=1}^{p+1}\frac{1}{p+1}\left(\genfrac{}{}{0ex}{}{p+1}{k}\right)B_{p+1-k}{n}^k\end{array}$

since only one term of the sum is modified. That is, a sum over one less integer can be written with the standard Bernoulli numbers without modification.

Bernoulli Polynomials

These polynomials are defined by the generating function

$\frac{t{e}^{xt}}{e^t-1}=\sum _{n=0}^{\infty }{B}_n\left(x\right)\frac{t^n}{n!}$

Inserting the definition of the Bernoulli numbers gives

$\sum _{n=0}^{\infty }{B}_n\left(x\right)\frac{t^n}{n!}=e^{xt}\sum _{k=0}^{\infty }{B}_k\frac{t^k}{k!}=\sum _{k,m=0}^{\infty }{B}_k{x}^m\frac{t^{k+m}}{k!m!}=\sum _{n=0}^{\infty }\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k{x}^{n-k}\frac{t^n}{n!}$

where the last rearrangement picks out individual powers of the variable t. The polynomials are thus given in terms of the numbers by

$B_n\left(x\right)=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k{x}^{n-k}=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_{n-k}{x}^k$

Again with the aid of a CAS or a bit of calculational persistence, the first few polynomials are

$\begin{array}{l}{B}_0\left(x\right)=1\\ B_1\left(x\right)=x-\frac{1}{2}\\ B_2\left(x\right)=x^2-x+\frac{1}{6}\\ B_3\left(x\right)=x^3-\frac{3}{2}{x}^2+\frac{1}{2}x\\ B_4\left(x\right)=x^4-2x^3+x^2-\frac{1}{30}\end{array}$

From the definition it is trivial to see that

$B_n\left(0\right)=B_n$

The value of the polynomials at unity is in most cases identical, which can be seen by shifting the quantity n in the all-important summation for Bernoulli numbers:

$B_n\left(1\right)=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k=\sum _{k=0}^{n-1}\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k+B_n=B_n$

This relation does not hold for $n=1$ because the intermediate sum is nonzero: instead one has $B_1\left(1\right)=-B_1$ .

The all-important summation leads to the following property of the polynomials for $n>0$ :

$\underset{0}{\overset{1}{\int }}dx{B}_n\left(x\right)=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k\frac{x^{n-k+1}}{n-k+1}{|}_0^1=\frac{1}{n+1}\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n+1}{k}\right)B_k=0$

The derivative of the polynomials is

$\frac{d}{dx}{B}_n\left(x\right)=\sum _{k=0}^{n-1}\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k(n-k)x^{n-k-1}=n\sum _{k=0}^{n-1}\left(\genfrac{}{}{0ex}{}{n-1}{k}\right)B_k{x}^{n-1-k}=n{B}_{n-1}\left(x\right)$

This relation holds for all values of the index and represents an Appell sequence. The relation can also be written in terms of an integral as

$\int dx{B}_n\left(x\right)=\frac{1}{n+1}{B}_{n+1}\left(x\right)$

with a shift in index. Constants of integration are here chosen proportional to Bernoulli numbers for consistency. Evaluating the left-hand side as a definite integral from zero to one is an alternate way to show that the polynomials have identical values at these two endpoints, apart from the exception noted above.

Euler-Maclaurin Summation

The formula derived here provides a way to evaluate the error between an integral and its approximation as a sum. While the derivation given in venerable Whittaker & Watson is accurate, it provides no hint as to why Bernoulli polynomials are employed as opposed to any other set of similarly defined polynomials. This derivation will be clearer in that crucial aspect.

An integral over a small region can be approximated to first order by a trapezoid,

$\underset{x_0}{\overset{x_1}{\int }}dxf\left(x\right)\approx \frac{x_1-x_0}{2}\left[f\right(x_0)+f(x_1\left)\right]$

where the heights at the endpoints are averaged and multiplied by the base. The error in this first approximation is clearly the difference of the two sides of this equation:

${\epsilon }_1=\frac{x_1-x_0}{2}\left[f\right(x_0)+f(x_1\left)\right]-\underset{x_0}{\overset{x_1}{\int }}dxf\left(x\right)$

Consider integration by parts of the following integral:

$\begin{array}{l}\underset{x_0}{\overset{x_1}{\int }}dx(x-c)f^{\prime }\left(x\right)=(x-c)f\left(x\right){|}_{x_0}^{x_1}-\underset{x_0}{\overset{x_1}{\int }}dxf\left(x\right)\\ =(x_1-c)f\left(x_1\right)-(x_0-c)f\left(x_0\right)-\underset{x_0}{\overset{x_1}{\int }}dxf\left(x\right)\end{array}$

Setting the constant equal to the average of the two endpoints, this is identical to the expression for error, which can thus be represented

${\epsilon }_1=\underset{x_0}{\overset{x_1}{\int }}dx(x-\frac{x_0+x_1}{2})f^{\prime }\left(x\right)$

Now make the change of variable

$x=(x_1-x_0)u+x_0=hu+x_0\frac{dx}{du}=h$

so that the integral becomes

${\epsilon }_1=h^2\underset{0}{\overset{1}{\int }}du(u-\frac{1}{2})f^{\prime }\left(x\right)$

where the argument of the function is left as is for conciseness.

The procedure going forward is to perform repeated integration by parts on this error integral to produce a sum of decreasing error terms. The polynomials that will appear during this process constitute an Appell sequence.

Most derivations at this point recognize the first factor in the integrand as the Bernoulli polynomial $B_1\left(u\right)$ and proceed without further comment. Such a procedure is problematic because this is also the Euler polynomial $E_1\left(u\right)$ , which differs from the Bernoulli polynomial by a constant of integration. Since constants of integration can be assigned arbitrarily in a definite integral, there needs to be justification for how they are chosen.

One way to do this is to focus on a final result that will be useful in evaluation. So far one has the exact statement

$\underset{x_0}{\overset{x_1}{\int }}dxf\left(x\right)=\frac{h}{2}\left[f\right(x_0)+f(x_1\left)\right]-h^2\underset{0}{\overset{1}{\int }}du(u-\frac{1}{2})f^{\prime }\left(x\right)$

over a small region. For extended regions one needs to add the contributions from each small region,

$\underset{a}{\overset{b}{\int }}dxf\left(x\right)=\frac{h}{2}\left[f\right(a)+f(b\left)\right]+h\sum _{k=1}^{n-1}f(a+kh)-\sum _{k=0}^{n-1}{h}^2\underset{0}{\overset{1}{\int }}du(u-\frac{1}{2})f^{\prime }\left(x_k\right)$

where $b=a+nh$ and $x_k=hu+a+kh$ . The first two terms represent the trapezoidal rule for estimating integrals: internal points contribute twice that of endpoints from adjacent small regions.

The third term also has contributions from adjacent regions, but given that it consists of integrals the contributions appear with opposite signs. If they can be made to cancel, the sum of error terms would then only depend on endpoint contributions. This will be the criterion for determining constants of integration.

Define an Appell-like sequence of polynomials with

$P_{i+1}\left(u\right)=\int du{P}_i\left(u\right)$

The integration by parts of a general error integral can then be written

$\begin{array}{l}{\epsilon }_i=h^{i+1}\underset{0}{\overset{1}{\int }}du{P}_i\left(u\right)f^{\left(i\right)}\left(x_k\right)\\ \phantom{{\epsilon }_i}=h^{i+1}[P_{i+1}\left(1\right)f^{\left(i\right)}\left(x_{k+1}\right)-P_{i+1}\left(0\right)f^{\left(i\right)}\left(x_k\right)]\\ \hfill -h^{i+2}\underset{0}{\overset{1}{\int }}du{P}_{i+1}\left(u\right)f^{(i+1)}\left(x_k\right)\end{array}$

Clearly the condition that allows cancellation of interior error contributions is

$P_i\left(0\right)=P_i\left(1\right)$

If this condition holds, the estimate to the integral becomes

$\begin{array}{l}\underset{a}{\overset{b}{\int }}dxf\left(x\right)=\frac{h}{2}\left[f\right(a)+f(b\left)\right]+h\sum _{k=1}^{n-1}f(a+kh)\\ -\sum _{m=2}^N(-1{)}^m{h}^m{P}_m\left(0\right)[f^{(m-1)}\left(b\right)-f^{(m-1)}\left(a\right)]+R\end{array}$

after repeated integration by parts. The remainder term R is simply the next error integral in the sequence and will not be considered further.

And now for the crucial part: determining the polynomials themselves. Starting from the one appearing in the integrand, the first few are

$\begin{array}{l}{P}_1\left(u\right)=u-\frac{1}{2}\\ P_2\left(u\right)=\frac{u^2}{2}-\frac{u}{2}+c_2\\ P_3\left(u\right)=\frac{u^3}{3!}-\frac{u^2}{2·2}+c_{2}u+c_3\\ P_4\left(u\right)=\frac{u^4}{4!}-\frac{u^3}{2·3!}+c_2\frac{u^2}{2}+c_{3}u+c_4\end{array}$

from which one can deduce the general form

$P_n\left(u\right)=\frac{u^n}{n!}-\frac{u^{n-1}}{2(n-1)!}+\sum _{k=2}^n{c}_k\frac{u^{n-k}}{(n-k)!}$

Equating the value of this polynomials and zero and one gives the equation

$\frac{1}{n!}-\frac{1}{2(n-1)!}+\sum _{k=2}^{n-1}\frac{c_k}{(n-k)!}=0$

With the choice $c_k=\frac{B_k}{k!}$ this becomes

$\frac{1}{n!}\sum _{k=0}^{n-1}\left(\genfrac{}{}{0ex}{}{n}{k}\right)B_k=0$

which is the all-important summation for Bernoulli numbers! Simply requiring these polynomials to have equal values at zero and one is enough to produce Bernoulli numbers and confirm that Bernoulli polynomials are indeed appropriate for the process. Comparing the general form here to the definition of Bernoulli polynomials gives

$P_n\left(u\right)=\frac{1}{n!}{B}_n\left(u\right)$

The approximation to the integral is thus

$\begin{array}{l}\underset{a}{\overset{b}{\int }}dxf\left(x\right)=\frac{h}{2}\left[f\right(a)+f(b\left)\right]+h\sum _{k=1}^{n-1}f(a+kh)\\ -\sum _{m=2}^N(-1{)}^m\frac{h^m{B}_m}{m!}[f^{(m-1)}\left(b\right)-f^{(m-1)}\left(a\right)]+R\end{array}$

and since the odd Bernoulli numbers in the sum are zero, this can be written a bit more simply as

$\begin{array}{l}\underset{a}{\overset{b}{\int }}dxf\left(x\right)=\frac{h}{2}\left[f\right(a)+f(b\left)\right]+h\sum _{k=1}^{n-1}f(a+kh)\\ -\sum _{k=1}^N\frac{h^{2k}{B}_{2k}}{\left(2k\right)!}[f^{(2k-1)}\left(b\right)-f^{(2k-1)}\left(a\right)]+R\end{array}$

The sum here is in general asymptotic, requiring an upper cutoff to prevent divergence.

It is instructive to compare this derivation to that originally given by MacLaurin as described in this paper, but with modern notation replacing fluxions. Start with the Taylor series of a function at some initial point:

$f\left(x\right)=\sum _{k=0}^{\infty }{f}^{\left(k\right)}\left(a\right)\frac{(x-a{)}^k}{k!}$

Integrating this from the initial point to one separated by unity will define a sum

$S_0=\underset{a}{\overset{a+1}{\int }}dxf\left(x\right)=\sum _{k=0}^{\infty }\frac{1}{(k+1)!}{f}^{\left(k\right)}\left(a\right)$

whose first term is $f\left(a\right)$ . Define a sequence of such sums with

$\begin{array}{l}{S}_1=f\left(x\right){|}_a^{a+1}=\sum _{k=1}^{\infty }\frac{1}{k!}{f}^{\left(k\right)}\left(a\right)\\ S_2=f\prime \left(x\right){|}_a^{a+1}=\sum _{k=2}^{\infty }\frac{1}{(k-1)!}{f}^{\left(k\right)}\left(a\right)\\ S_3=f^{″}\left(x\right){|}_a^{a+1}=\sum _{k=3}^{\infty }\frac{1}{(k-2)!}{f}^{\left(k\right)}\left(a\right)\end{array}$

where subscript on the sum indicates that its first term is $f^{\left(n\right)}\left(a\right)$ : initial constant terms cancel in the difference between endpoints, leaving only terms corresponding to the upper endpoint.

It is not hard to see from writing out a few terms in each sum that a particular derivative will appear only in a finite number of sums. Combining the sums with coefficients to be determined and picking out individual derivatives on the right-hand sides gives

$\sum _{m=0}^{\infty }{c}_m{S}_m=\sum _{k=0}^{\infty }{f}^{\left(k\right)}\left(a\right)\sum _{m=0}^k\frac{c_m}{(k+1-m)!}$

One can then solve for $f\left(a\right)$ by choosing the coefficients to make remaining sums on the right-hand side zero. If this is done, then one can write

$\begin{array}{l}f\left(a\right)=\frac{1}{c_0}\sum _{m=0}^{\infty }{c}_m{S}_m\\ \phantom{f\left(a\right)}=\underset{a}{\overset{a+1}{\int }}dxf\left(x\right)+\frac{c_1}{c_0}\left[f\right(a+1)-f(a\left)\right]\\ +\sum _{m=2}^{\infty }\frac{c_m}{c_0}[f^{(m-1)}(a+1)-f^{(m-1)}\left(a\right)]\end{array}$

Summing over n unit steps to a final point, one can write the integral as

$\begin{array}{l}\underset{a}{\overset{b}{\int }}dxf\left(x\right)=\sum _{k=0}^{n-1}f(a+k)-\frac{c_1}{c_0}\left[f\right(b)-f(a\left)\right]\\ -\sum _{m=2}^{\infty }\frac{c_m}{c_0}[f^{(m-1)}\left(b\right)-f^{(m-1)}\left(a\right)]\end{array}$

And how should the coefficients be chosen? With $c_m=\frac{B_m}{m!}$ the remaining sums are all

$\sum _{m=0}^k\frac{B_m}{m!(k+1-m)!}=\frac{1}{(k+1)!}\sum _{m=0}^k\left(\genfrac{}{}{0ex}{}{k+1}{m}\right)B_m=0$

which as anticipated is the all-important summation for Bernoulli numbers. The integral is thus

$\begin{array}{l}\underset{a}{\overset{b}{\int }}dxf\left(x\right)=\sum _{k=0}^{n-1}f(a+k)+\frac{1}{2}\left[f\right(b)-f(a\left)\right]\\ -\sum _{m=2}^{\infty }\frac{B_m}{m!}[f^{(m-1)}\left(b\right)-f^{(m-1)}\left(a\right)]\\ \phantom{\underset{a}{\overset{b}{\int }}dxf\left(x\right)}=\frac{1}{2}\left[f\right(b)+f(a\left)\right]+\sum _{k=1}^{n-1}f(a+k)\\ -\sum _{k=1}^{\infty }\frac{B_{2k}}{\left(2k\right)!}[f^{(2k-1)}\left(b\right)-f^{(2k-1)}\left(a\right)]\end{array}$

which as expected is equivalent to the result of the first derivation with $h=1$ .

While it may appear at first that Bernoulli numbers arise here in some sense from inverting a Taylor series, one can start the solution for derivatives at any arbitrary sum for the general result

$f^{\left(n\right)}\left(a\right)=\sum _{m=0}^{\infty }{c}_m{S}_{n+m}=\sum _{m=0}^{\infty }\frac{B_m}{m!}{S}_{n+m}$

Bernoulli numbers thus appear here, not from inverting a Taylor series per se, but from inverting its factorial coefficients. Curious...

Relations to Zeta Functions

Begin with the definition of the gamma function as an integral:

$\Gamma \left(z\right)=\underset{0}{\overset{\infty }{\int }}dt{t}^{z-1}{e}^{-t}$

Change the variable of integration using $t=(k+x)u$ ,

$\Gamma \left(z\right)=(k+x{)}^z\underset{0}{\overset{\infty }{\int }}du{u}^{z-1}{e}^{-(k+x)u}$

rearrange slightly and and sum over the introduced index,

$\sum _{k=0}^{\infty }\frac{1}{(k+x{)}^z}=\frac{1}{\Gamma \left(z\right)}\sum _{k=0}^{\infty }\underset{0}{\overset{\infty }{\int }}du{u}^{z-1}{e}^{-(k+x)u}=\frac{1}{\Gamma \left(z\right)}\underset{0}{\overset{\infty }{\int }}du\frac{u^{z-1}{e}^{-xu}}{1-e^{-u}}$

where the intermediate step uses an infinite geometric series. The left-hand sum is the Hurwitz zeta function $\zeta (z,x)$ expressed as a summation, and the right-hand side is its representation as an integral.

This integral converges as long as $z>1$ to avoid the singularity at the origin. For general values of the independent variable, interpret this integral as part of a complex contour that comes in from positive infinity, circles the origin and returns to positive infinity. That is, consider the contour integral

$\oint \frac{du}{u}\frac{(-u{)}^z{e}^{-xu}}{1-e^{-u}}$

The direct variable u has a complex argument of zero on the positive real axis. A contour that comes in from negative infinity has an argument of $-i\pi$ below the negative real axis, and an argument of $i\pi$ upon returning to negative infinity above the negative real axis. When this contour is rotated to handle the negative variable $-u$ these arguments are inverted, with the negative value on the part of the contour coming in from positive infinity above the positive real axis, and the positive value upon returning to positive infinity below the positive real axis.

Around the origin parametrize the contour as $-u=\rho e^{i\phi }$ . On this part of the contour the integral becomes

$\underset{0}{\overset{2\pi }{\int }}id\phi \frac{{\rho }^z{e}^{iz\phi }exp\left(x\rho e^{i\phi }\right)}{1-exp\left(\rho e^{i\phi }\right)}\approx {\rho }^{z-1}f\left(\phi \right),\rho \to 0$

where the left-hand side is exact, and approaches the right-hand value in the indicated limit. There is thus no contribution from this part of the contour, and the complete contour integral can be written

$\begin{array}{l}\oint \frac{du}{u}\frac{(-u{)}^z{e}^{-xu}}{1-e^{-u}}=\oint \frac{du}{u}\frac{(e^{\pm i\pi }u{)}^z{e}^{-xu}}{1-e^{-u}}\\ =e^{-i\pi z}\underset{\infty }{\overset{0}{\int }}\frac{du}{u}\frac{u^z{e}^{-xu}}{1-e^{-u}}+e^{i\pi z}\underset{0}{\overset{\infty }{\int }}\frac{du}{u}\frac{u^z{e}^{-xu}}{1-e^{-u}}\\ =2isin\pi z\underset{0}{\overset{\infty }{\int }}\frac{du}{u}\frac{u^z{e}^{-xu}}{1-e^{-u}}\end{array}$

where explicit values of the phase have been separated out in the second step as per the description above of the contour. Combining intermediate results gives

$\zeta (z,x)=\frac{1}{2isin\pi z\Gamma \left(z\right)}\oint \frac{du}{u}\frac{(-u{)}^z{e}^{-xu}}{1-e^{-u}}=\frac{\Gamma (1-z)}{2\pi i}\oint \frac{du}{u}\frac{(-u{)}^z{e}^{-xu}}{1-e^{-u}}$

using a reflection formula for the gamma function. One can now change the variable of integration using $u\to -u$ to get

$\zeta (z,x)=\frac{\Gamma (1-z)}{2\pi i}\oint du\frac{u^{z-1}{e}^{xu}}{1-e^u}$

where the contour is now rotated to the negative real axis, but still circles the origin. This representation is given here in a standard online reference.

Now insert the definition of Bernoulli polynomials,

$\zeta (z,x)=-\frac{\Gamma (1-z)}{2\pi i}\oint du\sum _{k=0}^{\infty }{B}_k\left(x\right)\frac{u^{k+z-2}}{k!}$

and let argument z be equal to a negative integer. Under the residue theorem, only terms with $k-n-1=0$ will contribute, giving

$\zeta (-n,x)=-\frac{\Gamma (n+1)}{2\pi i}\sum _{k=0}^{\infty }{B}_k\left(x\right)\oint du\frac{u^{k-n-2}}{k!}=-\frac{B_{n+1}\left(x\right)}{n+1}$

which is here in the standard online reference. Renaming the index gives the general relation

$B_n\left(x\right)=-n\zeta (1-n,x)$

Since the Hurwitz zeta function is analytic in both parameters, this relation can be used to extend Bernoulli polynomials as analytic functions of both index and argument.

When the second argument of the Hurwitz zeta function is unity, it becomes the Riemann zeta function. Since the value of Bernoulli polynomials at unity is usually the corresponding Bernoulli number, one has immediately

$B_n=-n\zeta (1-n),n>1$

This expression can appear with an additional alternating sign, which is unnecessary because the Bernoulli numbers that would be affected by such a sign are all zero.

Uploaded 2023.08.17 — Updated 2023.11.14 analyticphysics.com